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Abstract 

Nonlocal Hamiltonians are used widely in first-principles quantum calculations; the nonlocality stems 
from eliminating undesired degrees of freedom, e.g. core electrons. To date, attempts to couple nonlocal 
systems to external electromagnetic (EM) fields have been heuristic or limited to weak or long wavelength 
fields. Using Feynman path integrals, we derive an exact, closed-form coupling of arbitrary EM fields 
to nonlocal systems. Our results justify and clarify the couplings used to date and are essential for 
systematic computation of linear and especially nonlinear response. 

Theoretical study of the electronic response of quantum systems to electromagnetic (EM) probes is 
central to understanding experimental results and what they imply about the structure of matter. It is now 
possible to predict the response, from first principles, with sufficient accuracy to compare quantitatively 
to experiment E|. Nonlocal pseudopotentials are a key ingredient of many ab initio electronic structure 
calculations as they eliminate the largely inert and thus physically unimportant core electrons. The price 
paid is that the valence electrons feel a nonlocal atomic potential. 

Hence, it is vital to know how electrons moving in nonlocal potentials couple to EM fields. Early work 
found the coupling within the effective-mass framework. To date, refs. 0, [|, |j| represent the only first 
principles attempts to find the coupling, but with limited success. In ref. ||, the coupling terms involve 
an integral of the field over an unspecified path. Path ambiguity, discussed below, is a central problem 
that leads to non unique couplings and incorrect quantitative results. Ref. stops at the formal level due 
to operator ordering ambiguities, and an explicit, practical form of the coupling is given only for the long 
wavelength (i.e. dipole) approximation. This treatment excludes the description of AC magnetic fields or, 
more generally, spatially varying EM fields. Operator ordering problems also force the results of ref. || to 
be limited to linear (i.e. weak) coupling. 

Other attempts are also limited to linear coupling |?J where physical intuition leads to replacing the 
momentum by the velocity operator. When seeking nonlinear response, the coupling is either considered 
only for longitudinal fields ^ ^| or only within the long wavelength (dipole) approximation |jTof . 

Our work presents a rigorous derivation of an explicit, closed- form expression for the coupling of nonlocal 
systems to arbitrary EM fields. The value of this result is several fold: e.g., having the correct coupling- 
is essential for accurate and systematic ab initio calculation of linear and nonlinear response. Aside from 
first principles methods, many semi-empirical or phenomcnological models use nonlocal Hamiltonians with 
tight-binding like hopping terms. There have been only heuristic attempts to couple such systems to EM 
fields 0. Our results justify the assumed forms of previous couplings. 

Although ab initio calculations involve interacting electrons, nonlocal potentials are one-body operators. 
For clarity, we derive our results for a one-electron system. The results generalize directly to the interacting 
case. 

A main principle of electrodynamics is gauge-invariance. A particle in a local potential U(r) is governed 
by the Hamiltonian Hq(t,p) = p 2 /2m + U(r). The standard prescription for the gauge-invariant coupling 
uses the minimal substitution p — > II = p — A(f, t) (A is the vector potential). The coupled Hamiltonian is 
Ha(£, p) = II 2 /2m + U(r) + q(j)(r, t) (4> is the electrostatic potential). We use units where q/c = fi=\ (q is 
the particle's charge and c the speed of light). 
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The problem of interest, however, has a nonlocal potential V n i with Hamiltonian 

H = p 2 /2m + U(r ) + V nl , (1) 

where V n i(r, r') = (r|t^i|r') 7^ for r 7^ r'. Gauge-invariance demands that V n i depend on A. A major 
obstacle to finding H A is that V n i(r,r') is not specified in terms of the canonical operators (r, p). Below, 
we first rewrite V n i in terms of (f, p), and then perform minimal substitution in an unambiguous manner to 
find H A . 

Starting with V n i, we perform a Weyl transformation to obtain an equivalent operator W of (r, p) |T^ |: 

W(r, P) = Jds e"*-/ 2 V nl (r + |, r - |) e^' 2 . (2) 

It is straightforward to show that (r|VT(f, p)|r') = (r|V^j|r'). Thus V n i and W are the same operator. We 
now have written the Hamiltonian without fields as 

H Q (r,p) = p 2 /2m + U(r) + W(r,p). (3) 

For the equivalent classical system, W(r, p) is a function of the numbers r and p given by 

W(r,p) = /ds^(r + s/2,r-s/2)e^ ps . (4) 

At this stage, we would like to employ the minimal substitution p — > Tl in Eq. (^|) to find the coupled 
Hamiltonian H A . Unfortunately, the components of II do not commute so that initially equivalent orderings 
of operators in W lead to physically different couplings: e.g., p x p y = P y Px but Tl x tl y 7^ Tl y Il x since 
[n^jllj,] = iB 2 . Continuing along this line requires applying ordering postulates to the operators in W 
before minimal substitution to guarantee unambiguous final results. Only for long wavelengths (i.e. spatially 
uniform A) is the substitution unambiguous, and the results of |g Q are limited to this case. Similarly, 
this same ordering issue forced the results in S be limited to first order in A. 

An equivalent statement of the problem is that different evaluations of e~ l ^' s / 2 in Eq. (Q) correspond to 
translating by s/2 along different paths. This path ambiguity is present in ]3j where line integrals of A are 
used. When A is constant, all paths give identical results. But once A has spatial variation, different paths 
contribute different phases leading to non unique differing couplings. 

To find the solution, we first consider the classical system governed by W(r, p). Minimal substitution is 
well defined classically since p and r are commuting numbers. The classical coupled Hamiltonian is 

H A (r, p) = n 2 /2m + U(r) + W(r, II) + 0(r, t) . (5) 

Next, the Feynman path- integral formulation provides the quantum evolution of a system by summing over 
its classical dynamics. The propagation amplitude K between the two space-time points (r',0) and (r, t) in 
the operator formulation is given by Jf3| 

X(r,<;r',0) = (r|T 'exp (-% J dr H a (t)\ |r') 

= <5(r-r')-^<r|£ A (i)|r')+0(t 2 ). ( 6 ) 

Now K also equals the result of the Feynman formulation where we sum over all classical trajectories in 
phase space (r, p) weighted by each path's action S : 

K(r,i;r',0) = j' D[r] J D[p] e lS , (7) 

S = [ dr [p(r).r(r)-W A (r(r),p(r))] . (8) 
Jo 

The prime on the r path integral means that we only sum over paths where r(0) = r' and r(t) = r. 



2 



We seek the quantum Hamiltonian which is the 0(t) term in Eq. (^). Hence we will expand the 

path integral of Eq. (0) to 0(t) and extract £/a(*)- 

We now use the fact that W(r, p) is generally bounded, i.e. |W(r,p)| < C for some number C. This is 
true for atomic pseudopotentials V n i(r 7 r') which are smooth and have compact support. Thus when t — > 0, 
the contribution of the W term to S in Eq. (Q) is 0(t). Hence we need expand only to linear order in W: 

K = K loc + K nl , 



Kioc = J D[r] J D[p}e l 
K nl = -i J D[r] J D[p] J drW(r,n)e lS '° 



?/oc = 



J dr' [p • f - n 2 /2m - U(r) - 0(r, r')] ■ (9) 



The local propagator A"; oc describes the textbook situation of a local potential and is given by 

K loc = 6(r - r') - it (r\H loc \r') + 0{t 2 ) , (10) 

where Hi oc = Tl 2 /2m + U + q<f>. Thus, we concentrate on the nonlocal propagator K n i. 

The local potential U + q<fi depends only on r(r) and r and not on p. Hence as t —> 0, the contribution 
to Sioc from U + qcj) is 0(t). We ignore this contribution to K n i since it leads to terms of 0(t 2 ). However, 
p scales as p ~ (r — r')/t = 0{\/t) so expanding in t will fail for terms involving p. Thus we must perform 
the p integral explicitly. At fixed r, we change variables from p to II: 



K r , 



. J' D[r}e l f dT ' A ir E , (11) 
S = / ' D[U] J ^dTW(r(r),n(T))e l / dT ' [n -^ n2/2ml 
= JdrW (r(r), -^j J D[U y $ ^[n-r-tf /2 m] 

= N j ' drw[v(T),^-^e ml/2 ^ dT ' i2 . (12) 

Above, we use a functional derivative to pull W out of the II integral. This allows us to perform the Gaussian 
II integral, which yields a normalization constant M. 
Using the definition of W from Eq. (^) , we have 

i6\ /".,,/ s 



W {'-=W) ~ y*V.(r + |,r-|).-*. (13) 
The exponentiated functional derivative translates r: 

K nl = ~iN J dr Jds J D[v)Vm (r(r) + |,r(r) - |) X 

e i J dr'{ m [i- S S(T'-T)] 2 /2+A-i) ^ 14 ^ 



Note that changing the order of the integrals does not change the final result. We now change variables, at 
fixed s and r, from r(r) to f (r) as follows: 

r(r') = f (t') + [r' + s 9(t' - r) + (r - r' - s)r'/t] . (15) 



where d9{x)/dx = 8{x). This choice follows from the action of V n i in Eq. (14): V n i causes a transition 
from r(r) — s/2 to r(r) + s/2 at time r. We parameterize the paths as fluctuations f about a straight-line 
path with a discontinuity of size s at time r (term in brackets). With this choice, f obeys the conditions 
f (0) = f (t) = 0. 
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Upon substituting f , we have 

K nl = -iM Jdr j dsl 'r>[f]^„z(r(r) + |,r(r)-| 

g im/2(r-r'-s) 2 /*+i J dr' (mf 2 /2+ A-r) 

The conditions on f permit the Fourier expansion f(r') = ^ n f„ sin(7rnr'/t). The kinetic integral is 
J^dr't 2 = 7r 2 /2 J2 n n2 fn/t- We now note that for any analytic function g(x), J dx g(x)e tx ' > / t — y/iirt\g(0) + 
0(t 2 )\, so that within an integral we may set e™ 2 /' = ^/int[S(x) + 0(t 2 )]. Since K n i is already 0(£) and we 
ignore terms of 0(t 2 ), we replace e lx /* by \J int8{x) everywhere in the above integral: this applies to the 
f 2 /t as well as to the (r — r' — s) 2 /t terms. Physically, fluctuations in f and s contribute a wildly oscillating 
large phase as £ — > so only the stationary phase contribution at f = and s = r — r' remains when £ — > 0. 
Therefore, we set f = 0, s = r — r', absorb all constants from the integrations into Af' , and arrive at 

= -W [dTe*f dT ' A * x 



Km = 

Vm (x(r, r) + x(r, r) - ^ j , (16) 

x(r,r') = r(r')|f=o, 8 =r-r' = r' + (r - r') 0(r' - r) . 

Since we defined d9(x)/dx = 5(x) and S(— x) = S(x), we have 9(0) = 1/2 so x(r, r) = r' + (r — r')/2. 

The final hurdle involves evaluating J dr'A ■ x in Eq. (|l6|). We replace the <5(cc) appearing in x by any 
smooth, non- negative function g(x) of compact support about x = with width much smaller than £ and 
unit integral. We then replace 6(x) by G(x) where dG(x)/dx = g(x), so G increases monotonically from 
to 1 in a width much smaller than £. The integral is then 



J c£t'A(x(t,t'),t') -x(t,t') = ^ A(x,£) • dx + 0(£ 2 



We evaluate A at £ since any explicit time dependence of A changes JT n z by 0(£ 2 ), which is ignored. The 
integral is then just along a straight line from r' to r. 

After replacing the r integral by £ in Eq. (|l^), our nonlocal propagator is K ni = 
— itN 1 V n i(r, r') exp(j A • c£x) + 0(t 2 ). Recovery of the correct propagator for A = requires TV' = 1. 
Combining K n i with K\ oc (Eq. ([h])), we have as our central result the following gauge-invariant quantum 
Hamiltonian 

H A = (p-qA/c) 2 /2m + U + q4> + V£, (17) 
(r|t>„||r') = y ni (r,r')e^£>( x '*)- dx (straight-line). (18) 

In practice, we often use perturbation theory to compute response functions: we set Ha = Hq + H lnt and 
work to a desired order in the interaction H mt . We expand the above result to all orders in A: 



g(p ■ A + A ■ p) q 2 A 2 

2mc 2mc 2 



H A = H + q<j>- ^ n _ ^ + h~^ + V nT 



('l% n V> = Ka(r,r')X;i(^/ A-rfx) . (19) 



fc=i 



In relation to previous work, a coupling similar to Eq . (|l§| ) is given in Ref. Q but only after the straight 



line path is assumed with no justification. Refs. |J] and |1C both find a result equivalent to the special case 
of Eq. ( |l8| ) in the limit of spatially constant A. In Ref. the linear term (k = 1) of Eq. (19) is derived 
correctly, but attempts at finding a general expression for arbitrary A are stymied by operator ordering 
problems. 
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Figure 1: Absorption spectrum for diamond: e-iiui) vs. u> in eV. The dashed curve uses the momentum 
operator, the solid curve uses the velocity operator. 



We now apply our result to the case of linear optical response where the wavelength of A is much greater 
than atomic dimensions and we set A constant over the nonlocal range of V n i- For the Coulomb gauge 
V ■ A = (f> — 0, the linear coupling is (k = 1 term of Eq. (fl9|)) 

H{ nt = -A(i) -qv/c where v= -i[r,H ]/h. (20) 

Here, v = dx/dt is the velocity operator. This result justifies the intuition that A couples to the current 
j = gv. For a local system v = p/m as expected. The difference between velocity and momentum has been 
stressed when computing optical response |^, |^, [l0| [TBI]. However, those works assume the long wavelength 
limit or couple only to longitudinal fields. Our result is a direct demonstration that in the Coulomb gauge, 
the coupling of electrons to long waves is via the velocity operator. To show the significance of nonlocality, 
we perform ab initio pseudopotential density functional calculations |Hifl . Figure ^ shows the RPA absorption 
spectrum, £2(0;) with a 0.2 eV Lorentzian broadening, where either the momentum or velocity operator is 
employed. The effect of nonlocality is 15% to 20%, non negligible when aiming for quantitative comparison 
to experiment. 

Next, we study electrons in a uniform magnetic field. Magnetic fields provide a novel application of our 
result since A must vary spatially so that the usual long wavelength approximation is invalid. We calculate 
the magnetic susceptibility \ = —d 2 E/d 2 B of atomic carbon and neon. We choose A = — B(0,Q,x + a): 
the gauge center is offset from the ionic nucleus to simulate the realistic case when many atoms are present. 
We choose a = 4 a.u. as a reasonable value. We use s and p Kleinmann-Bylander nonlocal projectors Jlq] 
(p local) and compute \ following ref. 0. Calculating x requires the linear and quadratic couplings, both 
of which have nonlocal terms. For linear coupling, using v instead of p/m has been advocated 0, whereas 
the quadratic nonlocal term of Eq. ( [l9| ) is truly novel. Table ^ shows results for x when using various 
possible couplings: (1) the traditional local coupling (p — A) 2 /2m; (2) the linear coupling is from Eq. ( |l9| ) 
but the quadratic coupling is the local A 2 /2m; (3) both linear and quadratic couplings are from Eq. (|19|); 
(4) all-electron results excluding the contribution of the core Is states. Method (1) is qualitatively incorrect 
as it is not gauge- invariant. Using the correct linear coupling (method 2) greatly improves the quality of 
the results. However, the inclusion of our novel nonlocal quadratic terms (method 3) yields x m excellent 
quantitative agreement with the desired all-electron results. 

We emphasize that Eq. (|l8|) gives the coupling to any order in A. This allows, for the first time, for direct 
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Table 1: Negative of valence atomic magnetic suceptibility — x m cm 3 /mole. See text for details. 



and systematic calculation of high-order nonlinear responses without resort to longitudinal-only couplings 
H or to the long wave limit for the transverse case [p[ [l0[ p"7|. 

Finally, our result bears directly on systems with generic nonlocal hopping terms, e.g. tight-binding 
systems, where electrons hop from a localized state (3 on site R' to a state a at R with amplitude t a p(R, R'). 
Our result shows that, in the extreme tight-binding limit, the EM coupling modifies the amplitude to become 

i Q /3(R, R') exp Jj^ A • dx") . This justifies the straight-line integral used in the Peierls substitution and 



does not suffer from path ambiguity 1 1 1 
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